# Appendix J
# Santiago Lopez-Cariboni, 2022. "Political Regimes and Informal Social Insurance", Comparative Political Studies.

# ---- wd ----
wd <- 'replication_materials' 
setwd(wd)

# ---- packages ----
pkgs <- c("data.table", "foreign", 'panelView', 'fect', "lubridate", 
    "ggplot2", "dplyr", "mFilter", "texreg", "AER", "xts", "interflex", "gridExtra")
suppressPackageStartupMessages(sapply(pkgs,require,character.only=T))
# Plot the effects (use Johannes Karreth's function: https://github.com/jkarreth/JKmisc/blob/master/ggintfun.R)
# devtools::source_url("https://raw.githubusercontent.com/jkarreth/JKmisc/master/ggintfun.R")
source('code/functions/ggintfun.R')
# Function for time Lags
source('code/functions/lagpanel.R')


# ---- data ----
dt <- fread('data/dt_replication.csv')
dt <- dt[LDC==1, ]

# ---- Table 8 ----
dt <- dt[order(dt$iso3c, dt$year),] # sort for lagpanel
dt$l.r_cor <- lagpanel(dt$r_cor, dt$iso3c, dt$year, 1) 
dt$l.r_reg <- lagpanel(dt$r_reg, dt$iso3c, dt$year, 1) 



m <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * l.democracy
    + l.r_cor
    + l.r_reg
    - l.r_cor
    - l.r_reg
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt, na.action=na.omit)
 summary(m,  diagnostics = TRUE)
nc <- length(unique(m$model[['as.factor(iso3c)']]))

m.capacity <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * l.democracy
    + l.r_cor
    + l.r_reg
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt, na.action=na.omit)
summary(m.capacity,  diagnostics = TRUE)
nc.capacity <- length(unique(m.capacity$model[['as.factor(iso3c)']]))


m.capacity.controls <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * l.democracy
    + l.r_cor
    + l.r_reg
    + l.capacity
    + l.imports
    + l.exports
    + log(l.pop)
    + log(l.gdp.pcap)
    + l.elec.cons.pc
    + l.pop.density
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt, na.action=na.omit)
summary(m.capacity.controls,  diagnostics = TRUE)
nc.capacity.controls <- length(unique(m.capacity.controls$model[['as.factor(iso3c)']]))

v.names <- c(
    "Losses Cycle$_{t-1}$",
    "GDP Output Gap$_{t-1}$",
    "Democracy$_{t-1}$",
    "Democracy$_{t-1}$ $\\times$ GDP Output Gap$_{t-1}$",
    "Corporatization$_{t-1}$",
    "Independent regulatory agency$_{t-1}$",
    "State Capacity$_{t-1}$",
    "Imports$_{t-1}$",
    "Exports$_{t-1}$",
    "Population (log)$_{t-1}$",
    "Real GDP per capita (log)$_{t-1}$",
    "Electricity Consumption$_{t-1}$",
    "Population Density$_{t-1}$"
    )



table1 <- texreg(list(m, m.capacity, m.capacity.controls),
    file= "tables/OA/T8_capacity_elec.tex",
    label="tab:capacity_elec",
    custom.model.names = c('(1)','(2)', '(3)'),
    caption="Ciclicality of Electricity Losses in Autocracies and Democracies controlling for Electricity Sector Specific Measures of Enforcement Capacity",
    dcolumn = TRUE,
    no.margin=FALSE,
    fontsize="scriptsize",
    single.row=FALSE,
    use.packages=FALSE,
    booktabs = TRUE,
    digits=2,
    float.pos="htbp",
    sideways=FALSE,
    omit.coef= "(year)|(iso3c)|(Intercept)",
    include.rsquared = FALSE,
    stars = c(0.01, 0.05, 0.1),
    custom.gof.rows = list("Year dummies" = c("Yes", "Yes", "Yes"), 
        "Country fixed-effects" = c("Yes", "Yes", "Yes"),
        "Number of countries"= c(nc, nc.capacity, nc.capacity.controls)),
    # reorder.gof = c(1, 2, 5, 3, 4),
    custom.coef.names=v.names
    )




# ---- Table 9 ----
dt <- dt[order(dt$iso3c, dt$year),] # sort for lagpanel
dt$l.bti_mo <- lagpanel(dt$bti_mo, dt$iso3c, dt$year, 1) 
dt$l.milpercap <- lagpanel(dt$milpercap, dt$iso3c, dt$year, 1) 
dt$l.pubcorr <- lagpanel(dt$v2x_pubcorr, dt$iso3c, dt$year, 1) 


m <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * l.democracy
    + l.bti_mo
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt, na.action=na.omit)
summary(m,  diagnostics = TRUE)
nc <- length(unique(m$model[['as.factor(iso3c)']]))

m.c <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * l.democracy
    + l.imports
    + l.exports
    + log(l.pop)
    + log(l.gdp.pcap)
    + l.elec.cons.pc
    + l.pop.density
    + l.bti_mo
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt, na.action=na.omit)
 summary(m.c,  diagnostics = TRUE)
nc.c <- length(unique(m.c$model[['as.factor(iso3c)']]))


m5 <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * l.democracy
    + l.milpercap
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt, na.action=na.omit)
 summary(m5,  diagnostics = TRUE)
nc5 <- length(unique(m5$model[['as.factor(iso3c)']]))

m5.c <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * l.democracy
    + l.imports
    + l.exports
    + log(l.pop)
    + log(l.gdp.pcap)
    + l.elec.cons.pc
    + l.pop.density
    + l.milpercap
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt, na.action=na.omit)
summary(m5.c,  diagnostics = TRUE)
nc5.c <- length(unique(m5.c$model[['as.factor(iso3c)']]))


m6 <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * l.democracy
    + l.pubcorr
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt, na.action=na.omit)
 summary(m6,  diagnostics = TRUE)
nc6 <- length(unique(m5$model[['as.factor(iso3c)']]))

m6.c <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * l.democracy
    + l.imports
    + l.exports
    + log(l.pop)
    + log(l.gdp.pcap)
    + l.elec.cons.pc
    + l.pop.density
    + l.pubcorr
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt, na.action=na.omit)
summary(m6.c,  diagnostics = TRUE)
nc6.c <- length(unique(m6.c$model[['as.factor(iso3c)']]))


v.names <- c(
    "Losses Cycle$_{t-1}$",
    "GDP Output Gap$_{t-1}$",
    "Democracy$_{t-1}$",
    "Monopoly of Force (Bertlesmann TI)$_{t-1}$",
    "Democracy$_{t-1}$ $\\times$ GDP Output Gap$_{t-1}$",
    "(log) Military Personnel per 1,000 in pop (COW/WDI)$_{t-1}$",
    "Public sector corruption index$_{t-1}$",
    "Imports$_{t-1}$",
    "Exports$_{t-1}$",
    "Population (log)$_{t-1}$",
    "Real GDP per capita (log)$_{t-1}$",
    "Electricity Consumption$_{t-1}$",
    "Population Density$_{t-1}$"
    )


table1 <- texreg(list(m, m5, m6, m.c, m5.c, m6.c),
    file= "tables/OA/T9_capacity_disagregated.tex",
    label="tab:capacity_disagregated",
    custom.model.names = c('(1)','(2)', '(3)', '(4)', '(5)', '(6)'),
    caption="Ciclicality of Electricity Losses in Autocracies and Democracies controlling Dissagregated Measures of Capacity",
    dcolumn = TRUE,
    no.margin=FALSE,
    fontsize="scriptsize",
    single.row=FALSE,
    use.packages=FALSE,
    booktabs = TRUE,
    digits=2,
    float.pos="htbp",
    sideways=FALSE,
    omit.coef= "(year)|(iso3c)|(Intercept)",
    include.rsquared = FALSE,
    stars = c(0.01, 0.05, 0.1),
    custom.gof.rows = list("Year dummies" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), 
        "Country fixed-effects" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
        "Number of countries"= c(nc, nc5, nc6, nc.c, nc5.c, nc6.c)),
    reorder.gof = c(1, 2, 3, 5, 4),
    custom.coef.names=v.names
    )





# ---- Figure J1 ----
m_capacity_auto <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * l.capacity 
    + l.imports
    + l.exports
    + log(l.pop)
    + log(l.gdp.pcap)
    + l.elec.cons.pc
    + l.pop.density
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt[democracy==0,], na.action=na.omit)
pdf("figures/OA/J1_capa_auto.pdf", width = 8, height = 6)
ggintfun(obj = m_capacity_auto, varnames = c("l.outgap.gdp.hamilton", "l.capacity"), 
         varlabs = c("GDP Output Gap", "State Capacity"),
         title = FALSE, rug = TRUE,
         twoways = FALSE) +
coord_cartesian(ylim = c(-.15, .1), xlim =c(-1.8, 2.2)) +
ggtitle("Autocracies")
dev.off()


m_capacity_demo <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * l.capacity 
    + l.imports
    + l.exports
    + log(l.pop)
    + log(l.gdp.pcap)
    + l.elec.cons.pc
    + l.pop.density
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt[democracy==1,], na.action=na.omit)
pdf("figures/OA/J1_capa_demo.pdf", width = 8, height = 6)
ggintfun(obj = m_capacity_demo, varnames = c("l.outgap.gdp.hamilton", "l.capacity"), 
         varlabs = c("GDP Output Gap", "State Capacity"),
         title = FALSE, rug = TRUE,
         twoways = FALSE) +
coord_cartesian(ylim = c(-.15, .1), xlim =c(-1.8, 2.2)) +
ggtitle("Democracies")
dev.off()


m_bti_mo_auto <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * l.bti_mo 
    + l.imports
    + l.exports
    + log(l.pop)
    + log(l.gdp.pcap)
    + l.elec.cons.pc
    + l.pop.density
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt[democracy==0,], na.action=na.omit)
pdf("figures/OA/J1_monopoly_auto.pdf", width = 8, height = 6)
ggintfun(obj = m_bti_mo_auto, varnames = c("l.outgap.gdp.hamilton", "l.bti_mo"), 
         varlabs = c("GDP Output Gap", "Monopoly of the Force"),
         title = FALSE, rug = TRUE, rugsize = 0.2, jitter_factor = 2,
         twoways = FALSE) +
coord_cartesian(ylim = c(-.15, .2), xlim =c(0, 20)) +
ggtitle("Autocracies")
dev.off()

    
m_bti_mo_dem <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * l.bti_mo 
    + l.imports
    + l.exports
    + log(l.pop)
    + log(l.gdp.pcap)
    + l.elec.cons.pc
    + l.pop.density
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt[democracy==1,], na.action=na.omit)

pdf("figures/OA/J1_monopoly_demo.pdf", width = 8, height = 6)
ggintfun(obj = m_bti_mo_dem, varnames = c("l.outgap.gdp.hamilton", "l.bti_mo"), 
         varlabs = c("GDP Output Gap", "Monopoly of the Force"),
         title = FALSE, rug = TRUE,  rugsize = 0.2, jitter_factor = 2,
         twoways = FALSE) +
coord_cartesian(ylim = c(-.15, .2), xlim =c(0, 20)) +
ggtitle("Democracies")
dev.off()


m_r_prv_auto <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * r_prv 
    + l.imports
    + l.exports
    + log(l.pop)
    + log(l.gdp.pcap)
    + l.elec.cons.pc
    + l.pop.density
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt[democracy==0,], na.action=na.omit)

pdf("figures/OA/J1_militar_auto.pdf", width = 8, height = 6)
ggintfun(obj = m_r_prv_auto, varnames = c("l.outgap.gdp.hamilton", "r_prv"), 
         varlabs = c("GDP Output Gap", "Pivatization of State-Owned Utilities"),
         title = FALSE, rug = TRUE, rugsize = 0.2, jitter_factor = 3,
         twoways = FALSE) +
coord_cartesian(ylim = c(-.15, .1), xlim =c(1, 2)) +
ggtitle("Autocracies")
dev.off()


m_r_prv <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * r_prv 
    + l.imports
    + l.exports
    + log(l.pop)
    + log(l.gdp.pcap)
    + l.elec.cons.pc
    + l.pop.density
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt[democracy==1,], na.action=na.omit)
pdf("figures/OA/J1_militar_demo.pdf", width = 8, height = 6)
ggintfun(obj = m_r_prv, varnames = c("l.outgap.gdp.hamilton", "r_prv"), 
         varlabs = c("GDP Output Gap", "Pivatization of State-Owned Utilities"),
         title = FALSE, rug = TRUE, rugsize = 0.2, jitter_factor = 3,
         twoways = FALSE) +
coord_cartesian(ylim = c(-.15, .1), xlim =c(1, 2)) +
ggtitle("Democracies")
dev.off()



m_milpercap_auto <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * l.milpercap 
    + l.imports
    + l.exports
    + log(l.pop)
    + log(l.gdp.pcap)
    + l.elec.cons.pc
    + l.pop.density
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt[democracy==0,], na.action=na.omit)
pdf("figures/OA/J1_militar_auto.pdf", width = 8, height = 6)
ggintfun(obj = m_milpercap_auto, varnames = c("l.outgap.gdp.hamilton", "l.milpercap"), 
         varlabs = c("GDP Output Gap", "Military personnel per capita"),
         title = FALSE, rug = TRUE, rugsize = 0.2, jitter_factor = .1,
         twoways = FALSE) +
coord_cartesian(ylim = c(-.15, .1), xlim =c(0, 4.5)) +
ggtitle("Autocracies")
dev.off()


m_milpercap <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * l.milpercap 
    + l.imports
    + l.exports
    + log(l.pop)
    + log(l.gdp.pcap)
    + l.elec.cons.pc
    + l.pop.density
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt[democracy==1,], na.action=na.omit)
pdf("figures/OA/J1_militar_demo.pdf", width = 8, height = 6)
ggintfun(obj = m_milpercap, varnames = c("l.outgap.gdp.hamilton", "l.milpercap"), 
         varlabs = c("GDP Output Gap", "Military personnel per capita"),
         title = FALSE, rug = TRUE,  rugsize = 0.2, jitter_factor = .1,
         twoways = FALSE) +
coord_cartesian(ylim = c(-.15, .1), xlim =c(0, 4.5)) +
ggtitle("Democracies")
dev.off()


m_r_cor_auto <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * r_cor 
    + l.imports
    + l.exports
    + log(l.pop)
    + log(l.gdp.pcap)
    + l.elec.cons.pc
    + l.pop.density
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt[democracy==0,], na.action=na.omit)
pdf("figures/OA/J1_cor_auto.pdf", width = 8, height = 6)
ggintfun(obj = m_r_cor_auto, varnames = c("l.outgap.gdp.hamilton", "r_cor"), 
         varlabs = c("GDP Output Gap", "Corporatization of State-Owned Utilities"),
         title = FALSE, rug = TRUE, rugsize = 0.2, jitter_factor = 3,
         twoways = FALSE) +
coord_cartesian(ylim = c(-.15, .1), xlim =c(1, 2)) +
ggtitle("Autocracies")
dev.off()


m_r_cor <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * r_cor 
    + l.imports
    + l.exports
    + log(l.pop)
    + log(l.gdp.pcap)
    + l.elec.cons.pc
    + l.pop.density
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt[democracy==1,], na.action=na.omit)
pdf("figures/OA/J1_cor_demo.pdf", width = 8, height = 6)
ggintfun(obj = m_r_cor, varnames = c("l.outgap.gdp.hamilton", "r_cor"), 
         varlabs = c("GDP Output Gap", "Corporatization of State-Owned Utilities"),
         title = FALSE, rug = TRUE, rugsize = 0.2, jitter_factor = 3,
         twoways = FALSE) +
coord_cartesian(ylim = c(-.15, .1), xlim =c(1, 2)) +
ggtitle("Democracies")
dev.off()


m_r_cor_auto <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * r_reg 
    + l.imports
    + l.exports
    + log(l.pop)
    + log(l.gdp.pcap)
    + l.elec.cons.pc
    + l.pop.density
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt[democracy==0,], na.action=na.omit)
pdf("figures/OA/J1_reg_auto.pdf", width = 8, height = 6)
ggintfun(obj = m_r_cor_auto, varnames = c("l.outgap.gdp.hamilton", "r_reg"), 
         varlabs = c("GDP Output Gap", "Independent Regulatory Agency in the Electricity Sector"),
         title = FALSE, rug = TRUE, rugsize = 0.2, jitter_factor = 3,
         twoways = FALSE) +
coord_cartesian(ylim = c(-.15, .1), xlim =c(1, 2)) +
ggtitle("Autocracies")
dev.off()


m_r_cor <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * r_reg 
    + l.imports
    + l.exports
    + log(l.pop)
    + log(l.gdp.pcap)
    + l.elec.cons.pc
    + l.pop.density
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt[democracy==1,], na.action=na.omit)
pdf("figures/OA/J1_reg_demo.pdf", width = 8, height = 6)
ggintfun(obj = m_r_cor, varnames = c("l.outgap.gdp.hamilton", "r_reg"), 
         varlabs = c("GDP Output Gap", "Independent Regulatory Agency in the Electricity Sector"),
         title = FALSE, rug = TRUE, rugsize = 0.2, jitter_factor = 3,
         twoways = FALSE) +
coord_cartesian(ylim = c(-.15, .1), xlim =c(1, 2)) +
ggtitle("Democracies")
dev.off()






